Field theory for size- and charge asymmetric primitive model of 
ionic systems. Mean-field stability analysis and pretransitional 

effects 

A.Ciach 1 , W.T. Gozdz 1 and G.Stell 2 

1 Institute of Physical Chemistry, Polish Academy of Sciences 01-224 Warszawa, Poland 
2 Department of Chemistry, State University of 
New York Stony Brook, NY 11 794-3400, USA 
(Dated: February 6, 2008) 

Abstract 

The primitive model of ionic systems is investigated within a field-theoretic description for the 
whole range of size- , A, and charge, Z, ratios of the two ionic species. Two order parameters (OP) 
are identified, and their relations to physically relevant quantities are described for various values 
of A and Z. Instabilities of the disordered phase associated with the two OP's are determined in 
the mean-field approximation (MF). In MF a gas-liquid separation occurs for any Z and A / 1. In 
addition, an instability with respect to various types of periodic ordering of the two kinds of ions is 
found. Depending on A and Z, one or the other transition is metastable in different thermdynamic 
states. The instabilities found in MF represent weak ordering of the majority of the instantaneous 
states, and are identified with structural loci associated with pretransitional effects. 
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I. INTRODUCTION 



Universal features of collective phenomena can be determined within the framework of sta- 
tistical mechanics with the interaction potentials approximated by highly simplified generic 
models. The generic model that allows for a prediction of general features of the collective 
phenomena in ionic systems as diverse as molten salts, electrolytes, room-temperature ionic 
liquids, as well as systems containing charged nanoparticles (including proteins), and/or 
charged colloidal particles, is the primitive model (PM). In the PM ions are represented 
by charged hard spheres, and the solvent (if present) is taken into account only through 
the dielectric constant . In real systems the molecular structure of the solvent and 

multipole moments and polarizability of ions and solvent molecules affect the collective phe- 
nomena in the way that depends significantly on the specific properties of the ions and the 
solvent 0, [5], and is far from being understood. To identify and analyze specific effects 
in particular cases, however, one should know, in the first place, the phase diagram and 
structure in the generic case, where such effects are not present. 

Despite the simplicity of the interaction potential in the PM, the dependence of the phase 
diagrams and the correlation functions on the charge, 



and the diameter, 



A=£, (2) 



ratios of the two ionic species, respectively, is a largely unexplored problem. Until very 
recently the theoretical studies concentrated mainly on the restricted primitive model 

(RPM)Q,s,y,Q, where Z — X — I. In the RPM a gas-liquid separation occurs at low den- 
sities, and the bcc ionic crystal occurs at higher densities and sufficiently low temperatures. 
At still higher densities the fee crystal with and without substitutional order at low and at 
high temperatures respectively is stable 

Although much studied, the RPM is a quite artificial model, in the sense that no real 
ionic fluid has anions and cations of exactly the same size, although the size disparity is 
rather small in some, such as KCL. When one assumes full anion-cation symmetry, as in the 
RPM, there is a remarkable decoupling of the Ornstein-Zernike (OZ) integral equations that 
describe the density-density correlation function and the charge-charge correlation function, 
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respectively. The two correlation functions are still indirectly coupled through the "closure 
relations" that must be added to the OZ equations to yield a closed set of equations, and the 
different approximations that have come into standard use in solving the primitive-model OZ 
equations- the mean spherical approximation, the hypernetted chain approximation, etc- 
are defined by different closure relations [l[ • The total and the direct correlation function in 
the liquid theory are directly related to the (connected) correlation function and the vertex 
function respectively in the field theory (FT), and similar decoupling of the analogs of the 
OZ equations in the FT approach is found [3]. 

As soon as there is any degree of asymmetry, however, the density- density and charge- 
charge correlation functions are directly coupled, as are the OZ integral equations for these 



two functions as well as their analogs in the FT j^, 8]. We shall not go through the details 
of this case here, and the interested reader will find those details in Section 2.2 of Ref.jl]. 
There is one aspect of the results there that is worth noting here explicitly, however. If 
there is only a very small degree of asymmetry- for example, if the effective diameters of the 
anions and cations are very nearly the same in a 1-1 electrolyte, one would expect that the 
RPM equations would be "almost" satisfied in some sense, despite the direct coupling one 
finds in the asymmetric case. This turns out to be the case- one finds that deviations from 
the RPM result only appear very close to the critical point for very small asymmetries. The 



smaller the deviations, the c 
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24j results 



show that the PM with a small and moderate asymmetry exhibits qualitatively the same 



critical behavior as the RPM. In the asymmetric case the cp is s 
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23 



24]. 



and number densities (but higher volume fractions 
We expect to make further contact between our FT treatment for asymmetric ionic fluids 
and the general results of Ref.[l| concerning the critical behavior in future work. As far 
as we know, the effect of the size- and charge asymmetry on the gas-solid and liquid-solid 
transitions attracted much less attention so far. 



Despite our general insight into the asymmetric case that has come out of the work of 



Ref. 



and intensive studies of highly charged colloids 



25 



26. 



27 



28. 



29, 



30, £, 



32. 



33, 



34j . [35|, the case of the extreme asymmetry, Z, A — > oo, is not fully understood. This is 



because on the theoretical side the effect of approximations and assumptions on the results 



may be significant. In experiments 29 



30 



31 



36] the interactions differ from the 
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PM, and the role of specific interactions and/or of the solvent properties is not known in 
full detail. Moreover, results may be sensitive to impurities, nonequilibrium effects etc.. 
Last but not least, it is very difficult to obtain reliable simulation results for the phase 
diagram when Z, A —>■ oo . To ensure charge neutrality Z anions must be present per each 
cation, and in studies of collective phenomena the number of cations N + must be large. 
Typical experiments correspond to Z ~ 10 4 , i.e. to ~ 10 4 iV + ions in the simulation box. 
Because of the above reasons, there is no consensus concerning the phase behavior _m the 
PM in the case of the extreme asymmetry, and the results of some experiments 



are considered as controversial 



33 



30 



34] 



351 ] . The controversy concerns mainly the occurrence 



of the gas-liquid separation at room temperature in the case of monovalent counterions. 
Simulations indicate that at room temperatures the gas-liquid separation for monovalent 
counterions does not occur, but occurs for divalent counterions 33] . Preliminary field- 
theoretic results [8!] indicate that when the size and the charge of the cation are several 
orders of magnitude larger than the size and the charge of the anion, the gas-liquid type 
phase separation is preempted by a transition between gas and a colloidal crystal for a 
large range of temperatures, including room temperature p|. The predictions of Ref. [§] are 
limited to the mean field (MF) approximation, and the effect of fluctuations on the transition 
lines remains to be determined along the same lines as in the case of the RPM 37|. The 
phase behavior for higher volume fractions is also not quite clear. Many studies suggest 
a transition between the bcc and fee crystals, and next a gel or glass formation when the 



volume fraction of colloidal particles increases J3 



significantly influence the phase behavior 



ses [361 . 138 
ABB 



39] . Additional counterions and coions 
Note that the above phase 



behavior is quite different than the behavior in the RPM. 

The way the global phase diagram evolves when the asymmetry parameters are varied 
has not been systematically investigated yet. The values of Z and A have a strong effect 
on formation of ordered periodic structures of a crystal type. Coulomb interactions support 
structures in which the positive charges are compensated by the negative charges in regions 
as small as possible. On the other hand, precise compensation of charges in small regions 
may lead to an increase of entropy. The electrostatic and the entropic contributions to the 
grand potential both depend on the geometrical constraints. These constraints are due to 
a packing of Z spheres of a diameter cr_ and one sphere of a diameter cr+ = A<r_ into a 
charge-neutral aggregate, and a possibility of periodically repeating this aggregate in space. 
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Such an ordered structure is favored compared to the disordered phase when the entropic 
contribution is sufficiently small compared to the electrostatic energy. Combination of the 
geometric constrains and the charge-neutrality condition for some values of Z and A may 
promote ordering, while in other cases may suppress ordering. The competition between the 
periodic ordering and the phase separation also depends on Z and A. Thus, rich and complex 
phase behavior in the parameter space (£, T*; Z, A) may be expected, where £ and T* are the 
volume fraction of ions and temperature in the standard reduced (dimensionless) units 
defined in the next section. Depending on the values of Z and A, phase equilibria on the 
((,T*) phase diagram may include some (or all) of the following transitions: gas-crystal, 
gas-liquid, liquid-crystal, crystal-glass, crystal-gel or transition between different crystalline 
phases. For different Z and A the above phase equilibria may occur in quite different parts 
of the (C,T*) phase diagram. 

Some insight into the dependence of the phase behavior on Z and A can be gained 
from studies of pretransitional effects in the disordered phase. Such effects may be di- 
vided into two categories: effects associated with the gas- liquid separation and effects as- 



IE, 



15 



16 



41| emphasize the role 



sociated with the crystallization. Theoretical studies 

of cluster formation for the gas-liquid separation in ionic systems. Recent Monte Carlo 



(mc) da |20|, m 

dynamics (BD) 
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22,123 



24 



43 



44 



45| . molecular dynamics (MD) [46J and Brownian 



471 ] simulation studies provide interesting information on clustering in the 



size- and charge asymmetric PM. In Refs. 42 



43 



46 



471 | the simulations were performed at 



room temperature for several values of the volume fraction of the larger ions, for several sizes 
of both ionic species (size ratio 5-10), and for several values of the charge on the larger ion 
( 10-20 elementary charges). The cluster formation was studied for mono- di- and trivalent 
counterions for fixed values of the remaining parameters. In equilibrium the instantaneous 
structures are much more uniform in the case of monovalent counterions than in the other 
cases, where 'living' clusters containing two or more large ions are formed 



43 



47j. For 



trivalent counterions a substantial fraction of the clusters was found to be neutral. Attrac- 
tions between the neutral clusters are similar to the attractions between molecules, and may 
lead to the phase separation. The dependence on valency of counterions for a fixed charge 
on the larger ion and for fixed temperature is equivalent to a dependence on Z and T* . The 
question of the region in the parameter space ((,T*; Z, A) corresponding to a formation of 
neutral clusters remains open. 
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Pretransitional effects associated with weak long-range ordering are a subject of the 
present study. By analogy with the short-range ordering of the instantaneous states described 
above, we expect that such effects include long-range ordering of the instantaneous states. In 
the long-range ordered instantaneous states a local deviation from the uniform distribution 
of ions is periodically repeated in space. 

In this work we present an overview of the collective phenomena in the PM for the 
whole range of Z and A. Our FT analysis allows us to identify the dominant deviations 
from the uniform distribution of ions. We find structural lines that separate the stability 
region of the disordered phase on the ((,T*) phase diagram in the high-temperature part 
where the homogeneous instantaneous states are more probable than periodic states with 
small amplitude, and in the low-T* part, where the majority of the instantaneous states 
exhibit weak long-range order. At the structural line the amplitude of the periodic deviation 
from the uniform distribution of ions in the most probable instantaneous states increases 
continuously from zero. We stress that the structural line is neither a phase transition in the 
thermodynamic sense, nor a true spinodal. Determination of the probability of the periodic 
instantaneous states with large amplitudes, as well as determination of the phase transitions 
requires further studies, for which our analysis should serve as a starting point. 

In sec. II we describe the field theory for the PM. We introduce the concept of the struc- 
tural line and show that this line coincides with the MF approximation for the spinodal line. 
In sec. Ill general expressions for the eigenmodes and boundary of stability of the uniform 
phase in MF (structural lines) are given. The nature of the eigenmodes is discussed, and 
the dominant instantaneous states are identified in the same section. Sec. IV is devoted to 
the case of equal sizes. MF-spinodal (structural) lines associated with the gas-liquid sep- 
aration are described in detail in sec.V, and in sec. VI the structural lines associated with 
periodic ordering are discussed. Explicit results for representative values of Z and A are also 
presented in these two sections. A short summary is given in sec. VII. 
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II. FIELD-THEORETIC DESCRIPTION OF THE PRIMITIVE MODEL 



A. Coarse graining of the PM 

We consider a mixture of positively charged ions of a charge e + and diameter a + , and 
negatively charged ions of a charge e_ = — |e_| and diameter er_ in a structureless, incom- 
pressible solvent. In the PM the interaction potential of a pair a, (3, where the Greek indices 
denote + or — , is infinite for distances smaller than a a + <7g, i.e. we assume hard-core 
repulsions. The electrostatic potentials V a p{vi — r 2 ) between different pairs of ions a,j3 are 
the following 

V a0 (r) = ^O(r-a aP ), (3) 

where 

<T a fi = ^(<T a + <T{)) (4) 

and D is the dielectric constant of the solvent. The ^-functions above prevent from the 
contributions to the electrostatic energy coming from overlapping hard spheres, i.e. we do 
not include the electrostatic self-energy. 

In the field-theoretic approach we consider local densities of the ionic species, p+(r) and 
p_(r), and we include only smooth functions. For particular fields p+(r),p_(r) we assume 
that the electrostatic energy assumes the form 

U[p+,p-} = \j dr ±J dY 2 p a {r 1 )V aP (r 1 -r 2 )p /3 (r 2 ), (5) 

where summation convention for the Greek indices is used here and below. The prob- 
ability that the local densities assume a particular form p + (r),p_(r) is proportional to 
exp[— (3(U[p + , p_] — p+iV + — p_AL)], where N + = / drp + (r) and N_ = J drp_(r) are the 
number of positive charged and negative charged ions respectively. The chemical potentials 
of the two ionic species, p + and p_, are not independent, and have to be consistent with 
the requirement of charge neutrality, 

J drp + {v)e + = J drp_(r)|e_|. (6) 

Whereas the charge neutrality condition must be obeyed in macroscopic regions, it can be 
violated locally, in regions containing a small number of ions. The probability that the local 
densities assume the form p + (r),p_(r) is also proportional to the number of all microscopic 



states compatible with p + (r),p_(r). The number of all states can be written in the form 
exp(PTS), where (3 = (fcT) -1 , and where S, T and k are the entropy, temperature and 
the Boltzmann constant respectively. From the above it follows that the Boltzmann factor 
assumes the form 

p[p + ,p_] = S- 1 exp(-^[p + ,p_]), (7) 

where 

~=Jd P+ J Dp_exp(-(3n MF [p + ,p.}) } (8) 

and 

n MF [p + ,p_] = F,[p + ,p_] + f/[p + ,p_] -p + J rfrp + (r) -p_|drp_(r) (9) 

is the grand potential in the system where the local concentrations of the two ionic species 
are constrained to have the forms p + (r),p_(r). In Eq.([9]) Fh = — TS is the Helmholtz free 
energy of the hard-core reference system. In the local density approximation i 7 ) l [p + ,p_] = 
/ drfh(p + (r), p_(r)), and fh consists of the ideal-gas contribution plus the excess free-energy 
density of hard-spheres with different diameters f^ x , 

/3A(p + (r),p_(r)) = £p Q (r)(ln(A3p Q ( r )) - 1) + (p + (r), p_(r)). (10) 

a 

The A a is the thermal de Broglie wavelength of species a, a = +, — . Here for f^ x we 



assume the Percus-Yevick compressibi 



ity route to the Helmholtz free-energy of a hard- 



sphere mixture with size asymmetry 48(]. The above theory is equivalent to an approximate 
version 49] of the exact collective- variables theory [50I . 51]. 

As a length unit we choose The dimensionless densities and inverse temperature 

are chosen as 



4 = Af*+-. f-=P-o%-, P* = P^^- (11) 

DO _| 

As thermodynamic variables we choose T* = 1//3* and the volume fraction of all ions, 



C = ^(p +4 + Po-^). (12) 

The volume fraction of all ions is related to the quantity proportional to the number density 
of all ionic species, 

s = ^(Po+ + Pi-) = «i-(JV+ + N-)/V, (13) 
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where V is the volume and N± denote the average numbers of positive and negative charges 
for given chemical potentials. The relation between ( and s is given by 

t = s[l + 36{6-u)-6 3 v]. (14) 

We also define the parameters describing the charge- and size asymmetry, 

z-i 

v = —y < 15 > 

and 

8 = ^. (16) 

A + 1 v ; 

For 8v > (8v < 0) the charge at the larger ion is larger (smaller). Because of the symmetry, 
we limit ourselves to 5 > and consider —1 < v < 1. 

Let us focus on homogeneous instantaneous states corresponding to the extremum of 
Q MF . The densities in the uniform states, p^,(r) = const., corresponding to 5VL MF / 5p* a = 0, 
will be denoted by Pq q . At sufficiently high T*, homogeneous instantaneous states correspond 
to the global maximum of the Boltzmann factor ([7]), i.e. to the global minimum of Q MF . 
At low temperatures the extremum of Q MF may correspond to either a local minimum, 
a saddle point or to a maximum. For small deviations of p+(x) and pl(x) from and 
Po- respectively, the functional (02) can be expanded. The expansion can be truncated for 
P* a — Poa ~^ 0- We write Q MF in the form 

n MF = n + n 2 + n inU (17) 

where Q = Q mf [pq + , p$_] is the value of the functional at the extremum, and Q 2 denotes 
the Gaussian part of the functional. In Fourier representation we have 

1 r r/k 

/?fi 2 [Ap + (k), Ap_(k)] = -] j^Ap a (-k)C° a ,(k)Ap,(k), (18) 

where Ap a (k) is the Fourier transform of Ap* (r) = p* (r) — pl a . In Eq.(fl8l) and hereafter 
the wave numbers are in a+i units. The second functional derivative of Q MF [p+,p*_] at 
Pa = Poa consists of two terms, 

C° a Jk) = a aP + (3V aP (k). (19) 



The first term is given by 



a " = 9^,' <20) 
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where the derivative is taken at p* = p^ a . Explicit expressions for a a p are given in Appendix 
A. We have used the chemical potentials for the asymmetric hard-sphere mixture obtained 
in the Percus-Yevick approximation in Ref. 48j. The second term in (7?g(k) is the Fourier 
transform of the potential ([3]) and we find 

r 47i cos(kr + ) 



(3V ++ (k) = (3*Z 
/3VL_(Jfe) = (3*Z 



k 2 

1 4tt cos(Ar_ 



PV+-(h) — —(3 : 



k 2 
47T cos k 



where 



<?+- 



k 2 



1±5. 



The higher-order part of the functional (IT7I) can be written in the form 



where 



dr 



"3l 



Ap Q (r)A P/3 (r)Ap 7 (r) + ^Ap a (r)Ap,(r)Ap 7 (r)A Pi ,(r) 



4! 



da, 



a/3 



and 



da, 



dp* af3lu dp* v 

and the derivatives are taken at p* = Pq q . 

The RPM limit of the above general model corresponds to 8 — v 



a/37 



(21) 

(22) 
(23) 

(24) 

+ ..., (25) 
(26) 

and was studied 



before in Ref. 



The opposite, col 



loid limit of 5 -> 1, v -> 1 (A, Z -> oo) with s = 0{Z°) 



(i.e. ( = 0{Z x )) is described in Ref. 



B. Correlation and vertex functions and their generating functionals 

In the coarse-grained description the densities p*(r) play the role of the microstates. 
Thus, the standard definitions of thermodynamic potentials and correlations functions can 
be applied. In particular, the grand-thermodynamic potential in a presence of external fields 
J + (r) and J_(r) is given by 

fi[J+(r),J_(r)] = -mogS[J + (r),J_(r)]. (27) 

where 

S[J+(r),J_(r)] = ^ -Dp+ | J Dp-exp[-/3(fi MF [p + ,p_]-yj Q (r)p Q (r))]. (28) 
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The external fields may play an important role in some experimental cases. For systems 
where no external fields are present, as in our case, they are introduced as auxiliary fields 
for computational reasons. The Legeandre transform 

F[p + (r),p_(r)] = Q[J + (r), J_(r)] + jf J«(r)p a (r), (29) 

where p a = — 6Q[J+(r), J_(r)]/5J a , is the density functional. The — f2[J + (r), J_(r)] and 
F[p + (r), p_(r)l are the generating functionals for the correlation and the vertex functions 
respectively [53j, |54|. The two-point vertex function C Q/ 3 = 5 2 F/5p a 5p l 3 is related to the 



analog of the direct correlation function c a p. In Fourier representation this relation has the 
form [531] 

C a/3 (k) = f^-c a p(k), (30) 

where d^s is the Kronecker symbol, and (p a ) is the average density of the specie a. In terms 
of the vertex functions the FT analogs of the OZ equations take the simple form 

C a p{k)G 01 {k) = (31) 

where the Fourier transform G a p(k) of the two-point correlation function, 

s 2 n 



G a0 (r) = <Ap Q (0)A P/3 (r)> 



5J a (0)5Mr) 



is related to the analog of the total correlation function h a p by 53] 



(32) 



G a p(k) = 6$(p a ) + h aP (k)(p a )(pp). (33) 

The FT approach with the local-density approximation (flOj) is designed for a description of 
the long range ordering that determines phase transitions. Due to the coarse-graining (or 
'smearing' of the hard-spheres structure) the correlation functions for distances < a + + a_ 
are not correctly reproduced. 



C. Mean-field approximation 

In practice the exact forms of Q and F cannot be obtained. In order to adopt standard 
approximations to the PM, let us focus on the most probable instantaneous distributions of 
ions Poa( r )> corresponding to the global maximum of the Boltzmann factor (in the presence 
of the external fields the latter is proportional to the integrand in Eq. fl28l) ). At low T* 
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the most probable distributions may not be uniform; in particular, crystalline phases are 
characterized by distributions of ions that are periodic in space. We rewrite Eqs.( l2"7|) and 
in the form 



n[J+, J-} = n MF [p* 0+ (r),p*_(r)} - J Ja(r) P * 0a (r) - kT\ogA~[J +1 J 



(34) 



and 



F[p+,p_] = n MF [p* 0+ (v),p* _(v)] + jf J a (r)(p„(r) - p* 0a (r)) - kT log AE[J + , J_], (35) 
where 



AH[J+, J_] = y DAp+ J DAp_ exp [ - (3[An MF [A P+1 Ap_] - jf J a (r) Ap 



(36) 



and 



AO*" [Ap + , Ap_] = ^ [p* + + Ap + , p*_ + Ap_] - [p* + , p*_]. (37) 

The last term in Eq.( l34l) is the contribution to the grand potential associated with fluc- 
tuations Ap + , Ap_ around the most probable distributions, Po+( r ) an< ^ Po-( r )- m the MF 
approximation this term is just neglected, and the grand thermodynamic potential is ap- 
proximated by the minimum of fl MF [p + , p_] — f r J a (r)p a (r) with respect to p a for fixed 
J a (r). The average values of the local densities, 



(p Q (r)> = p a (r) = p* Q (r) + A~ [J + , J^ 1 J DAp + J DAp_Ap a {v) x 

p(^AQ MF [Ap + , Ap_] - f J Q (r)Ap Q (r 



(38) 



cxp 



in MF are identified with their most probable values Po a (r), i-e. the fluctuation contribution 
in Eq.(l38l) is neglected. Bey ond MF (p Q (r)) and Po a (r) may differ from each other, as found 



in the case of the RPM 



37 



551 ] . In the context of the PM the role of the fluctuation contri- 



bution in 
approach 



5qs. and fl38l will be studied in the future work, by following the Brazovskii 
5J| applied to the RPM in Ref. 



4l|. 



In the MF approximation the generating functionals for the vertex and the correla- 
tion functions, F[p + ,p_] and — Q[J + ,J_], reduce to f2 MF [p + ,p_] and — (jl MF [p + , p_] — 
J T J a (r)p a (r)^j respectively. For J a = the two-point correlation and vertex functions, 
G° a p and C£g respectively, in the above version of MF are identical to the corresponding 
functions in the Gaussian approximation. In the Gaussian approximation Q MF in the Boltz- 
mann factor (jTJ) is approximated by Eq. (1171) with the last term neglected. For Q MF = Q 2 
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the exact two-point correlation functions, G^p, are related to given in Eq. (fT9|) by the 
OZ equations analogous to Eq. flBTj) . 



D. Phase transitions, spinodals and pretransitional effects 

Let us focus on the FT predictions for the phase diagram in the case of J a = 0. The 
stable (metastable) phases correspond to the global (local) minimum of F[p +1 pJ\ (Eqs. fpZUj) 
and (1351) ). Two minima corresponding to different forms of p a (r) are of equal depth at the 
coexistence between the corresponding phases. Let us assume that the global minimum of 
F[p + , p_] corresponds to Po a (r), which deviates significantly from the constant function p^ a , 
corresponding to another extremum of F[p +1 pJ\. Crystalline phases are characterized by 
such density distributions. The minmum of F corresponding to such nonuniform distribu- 
tions can be found only when Q int in Eqs.(TT71) and (131)1) is included, because the Gaussian 
part 1^2 describes correctly only infinitesimal deviations from the uniform distributions of 
ions Poa( r )- ^int is given in terms of the three- and higher-order vertex functions in the MF 
approximation. Thus, the two-point vertex (or correlation) functions are not sufficient to 
find the transitions to the ordered phases; it is essential to include in Eq.fl35l) the higher-order 
vertex functions as well. 

Let us consider the boundary of the global (local) stability of the uniform phase. The 
extremum of F changes character at the corresponding continuous-transition (spinodal) 
line. At the high and at the low-T* side of this line the homogeneous phase corresponds 
to a global (local) minimum, and to a saddle point of F respectively. More precisely, for 
given C, the boundary of stability of the uniform phase is given by the highest T*, for which 
det C a p(k) = (5F/5p a 5p/3 is not positive definite) for some k = In MF the boundary of 
the stability of the unifo rm p hase is obtained with C a p(k) approximated by C®Jk) (Eq. (IT^l) ). 



From the earlier studies 



41 



56| it follows that when det C%(k) = for k ^ 0, the fluctuation 
contribution to F in Eq. fl35|) is comparable to Q MF . Hence, the difference between C a p(k) 
and C®p(k) may be comparable to C^Jk). In such a case the fluctuation contributions to 
Eqs. fl35|) - fl38|) lead to significant shifts of the phase boundaries. In addition, if a continuous 
transition to a phase with periodic densities is found in MF, a fluctuation-induced first- 



order transition is expected at significantly lower temperatures 
contribution is included. 



41 



56l | when the fluctuation 
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The analysis of the stability of Q MF [p + , p_] cannot give information on the actual in- 
stability of the disordered phase, but may give information on the pretransitional effects 
associated with weak long-range ordering of the instantaneous states. Let us compare the 
probability that an instantaneous state Ap + (x) 7^ 0, Ap_(x) 7^ occurs, with the probabil- 
ity of finding the homogeneous instantaneous state Ap + = Ap_ = 0. According to Eqs.([7J) 
and (TTT|) . the probability ratio is 

p[Ap + , Ap_] = p {Q2+Q . nt) 

p[0,0] ' 1 ; 

Let the smaller eigenvalue of C^Jk) in Eq. (fl8l) be denoted by C^Jk). The eigenmode 0(k) 
corresponding to k = ko can be written as a linear combination of the functions 

&(k|*b) = ^t(w5{\t - k 0t ) + w*6(k + k oi )) (40) 

with different directions of the vectors koi, where |koi| = &o and ww* = 1. Let us focus on 
the instantaneous structure that has the form 0(k) = $^j(k|A;o), and the amplitude of the 
other eigenmode of C^Jko), ^(k), vanishes. From Eq. (fT8l) we have (3VL 2 = VJ2i ^fC^ka) /2, 
where V is the volume of the system. Then Eq. (l39|) takes the form 

Ti $ 2 V / ~ 



p[$^(k|A; ),0] 

- = exp 



C H (k ) + O(^)) 



(41) 



p[0,0] 

In the above p is the probability (1391) as a function of the new variables, 0(k) and fj(k). 

The sign of the term 0($j) depends on the form of VLi nt in Eq. (II 71) . and on the form 
of </>(k). Here we limit ourselves to infinitesimal amplitudes When the term 0($j) in 
( )4T|) is negligible, the probability ratio is larger or smaller from unity for C^(fco) < or for 
C<f>(l>{ko) > respectively. For given (, the solution of C^ko) = yields the temperature 
^*(C> ^0); an d C^ko) > is equivalent to T* > T*((; k ) (see ( fl9l) and below). The most 
probable eigenmode corresponds to k b such that C^{k) assumes a minimum at k — kh. 
Hence, the instantaneous uniform state occurs with a higher probability than the periodic 
instantaneous state with infinitesimal amplitude and any wavelength for C^kf,) > 0, i.e. 
for T* > T*(C; kb). When for given ( temperature decreases from T* = T*(£; kb), the range 
of k such that p[$j^i(k|/c), 0]/p[0, 0] > 1 increases. This means that the population of the 
instantaneous periodic states that occur more frequently than the instantaneous uniform 
state increases. Ordering of the significant fraction of the instantaneous states can be recog- 
nized as a pretransitional effect associated with the transition to the corresponding ordered 
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phase. The line in the phase diagram corresponding to C^(A;&) = can be identified with 
the structural line, because different instantaneous states dominate on different sides of this 
line. Note that the structural line coincides precisely with the boundary of stability of the 
uniform phase in the MF approximation (MF spinodal). 

Near the minimum of C®p(k) the probability ratio in Eq. (j4ip varies slowly with k. Further 
studies are required to determine the probability that different wave packages corresponding 
to k ~ fcfe are thermally excited, and to verify if the wave packages are related to clusters of 
ions in real space. 

We should note that the large-amplitude periodic states may be more probable than the 
uniform state even for C®p{kb) > 0, if the term in Eq. (1411) is negative. Note also that 

in the above studies we assumed that the instantaneous states did not include the eigenmode 
r) associated with the larger eigenvalue of C^Jk). Interesting information on pretransitional 
effects can be gained from studies of the effect of coupling between the two eigenmodes 
in Clint- In particular, the probabilites of the eigenmodes rj, and 0, given by J Dct)f}[4>, v] 
and / Drjp[(j), rj\ respectively, are of considerable interest. In the case of the RPM [57] the 
charge- charge correlations play a role analogous to short-range interactions between the 



ions 



57|, 



581 ]. and at low volume fractions / D<pp[<f>,T]] turnes out to be important for the 



phase separation [41|, |57| • The above questions go beyond the scope of this work. 



III. EIGENMODES AND MF BOUNDARY OF STABILITY OF THE DISOR- 
DERED PHASE IN THE GENERAL CASE OF ARBITRARY 5 AND v 

A. Eigenmodes 

Our purpose here is to find the instantaneous states that may lead to instabilities of the 
disordered phase for arbitrary 5 and v in the MF. In order to find the spinodal lines in 
the MF (structural lines), it is sufficient to analyze /3f2 2 (see Eqs.(IT8l) and (1391)). C^Jk) 
(Eq. (fl9l) ) can be diagonalized, and S7 2 can be written in the form 

^ = 2 / (2^ K^W'K-k) + mC° m (k)v(-k)} . (42) 
Both C^Ak) and C~Jk) depend on the wavenumber k in a nontrivial way, 



d o {k) = C\ + {k) + C°_(k) ~ szgn(C° + _(k))B(k) 
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and 

d o {k) = C° ++ (k) + C°-(k) + sign{C°_(k))B(k) ^ {u) 

where 



B(k) = \]A 2 + 4C£_(£;) 2 , (45) 

and 

A(k) = sign{Cl_{k))[C°__{k) - C° ++ (k)]. (46) 

In the above C®p are given in Eq. (115)1 . a a p are given in Appendix A and V a p{k) are given in 
Eqs.(|2IJ-P3J. 

In the general case the corresponding eigenmodes have the forms 

0(k) = 5(*)Afl.(k) - &(*)A£. (k), (47) 
7j(k) = 6(fc)Ap;(k) + a(A;)Ap* (k), (48) 



with 

g(Jfe) _->K*) + B(*)^ 



b(k) 



2B{k) 



-A{k) + B(k) 1/2 



(49) 



(50) 



2B(k) 

The eigenmodes represent two order-parameter (OP) fields, and in principle either one of 
them may lead to instability of the disordered phase. 

On the MF level the FT analogs (1ST)) of the OZ equations are decoupled in the eigenbasis 
of the vertex functions, i.e. 

G° H (k) = l/C° M (k) and G° vv (k) = 1/C° m (k). (51) 

Note that in the fully symmetrical case {y — 6 — 0), C°__(k) = C+ + (k), and a(k) = b(k) = 
1 / \/2, hence and r] are proportional to the charge- and the number density respectively 
Moreover, Eqs.gSJ and flH} reduce to C° H (k) = (Cl + (k) + C°_(k) - 2C°_(Jfe))/2 and 
Crm(k) = (C++(k) + C°__(k) + 2C ( l_(k))/2, representing the charge-charge and the density- 
density vertex function respectively. When any asymmetry is present and C°__(k) ^ C+ + (k), 
CL, and C® v that satisfy the decoupled OZ equations differ from the charge- and number 
density vertex functions respectively, in agreement with Ref.[l|. 
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B. Dominant order parameters 



Let us focus on thermodynamic conditions such that the uniform distributions p$ a corre- 
spond to the minimum of Q MF , and consider the two OP fields with equal wavelengths and 
equal infinitesimal amplitudes. The field leading to a smaller increase of Q MF is thermally 
excited with a higher probability ([7]), and such OP dominates over the other one. Which one 
of the two infinitesimal fields, 0(k) or r/(k) dominates, depends on which function, C^k) 
or is smaller (see (H2j) ). This, in turn, depends on sign{C + -{k)) (see (|4~3]) and (pf4l) ). 

From Eqs. (fT9]) and (123]) we obtain the line 

mit 4tt cos k i . 
T* = —p-a£ (52) 

separating the phase space ((, T*) into the high-temperature part where the OP 0(k) dom- 
inates, and the low-temperature part where the other eigenmode dominates. In the above 
is a function of Q (or s) given in Appendix A. Note that for different wavenumbers k 
Eq. (l52"I) yields quite different lines. In particular, for long-wavelengths k — > 0, correspond- 
ing to phase separation into two uniform phases, from ( l52l) we find that ?7(0) dominates 
for all temperatures T* < oo. We can thus conclude that if the phase separation (stable or 
metastable) occurs, it is induced by the field rj. On the other hand, for k > it/2, correspond- 
ing to the periodic ordering of ions with the wavelength 27r/k < 4 (in units), the OP 

0(k) dominates for all temperatures T* > 0. In principle, the uniform system may become 
unstable with respect to the eigenmodes with 0<A;<7r/2as well. As we show later, this 
is indeed the case for very large asymmetry. In such a system 0(k) or fj(k) dominates for 
temperatures higher or lower than the temperature given in Eq. (l5"2l respectively. 

C. The MF spinodal lines 

Let us first analyze the stability of the uniform phase with respect to long- wavelength 
deviations from homogeneous distributions of ions in the MF approximation. For k — > 
Eqs.f l44p and (T4"3"|) assume the asymptotic forms 

cL(k) = a ++ z " + a -Z + ^ + --^ 2 + 0{k% (53) 



and 



C^W- ^^'W (54) 
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C^ik) diverges for k — > 0. When 0^(0) = oo, then the probability p oc exp(— (3Q MF [(j), rf\) 
of thermally exciting any nonzero 0(0) vanishes (see (IT7|) and fj42|) ). as required by the global 
charge neutrality (7° (A;) can vanish for k = when 

The above equation describes the MF spinodal line associated with the phase separation 
into two uniform phases characterized by different values of the OP r\. 
The MF instability with respect to the fc-mode is given by 

det<^(*) = <^(A0<^(A0 = O, (56) 

where 

C° m (k)C%(k) = (a ++ + (3V ++ (k))(a^ + pV—{k)) - (a + _ + (3V + 4k)) 2 = (57) 

sm 2 (k5). A _. 2 b(k, s, 5, v) nit 



k 4 K r ' k 2 

and the functions b(k, s, 5, v) and d(s, 5, v) > are given in Appendix B. 

Note that there is a single positive solution of Eq. (|56|) for (3* for all values of the remaining 
parameters. The wavector at the MF spinodal (structural line) is determined from 

d(C° m (k)C° H (k))/dk = 0. (58) 

The solutions of Eqs. fl56|) and fl58|) give the spinodal line and the wavenumber of the 
dominant instantaneous state. 

D. Nature of the eigenmodes for different 5 and v 

It is instructive to discuss the nature of the eigenmodes that may drive the system out 
of the disordered state. The equations (I47p and (1481) can be easily inverted to give 

Apl(k) = 5(A;)0(k) + 6(A0?ftk), (59) 

Ap* (k) = -b(k)4>(k) + d{k)fj(k). (60) 

The local charge density, 

q(r) = ZAp* + (r)-Ap*_(r) (61) 
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(in |e_|<7 + ^ units), and the local deviation of the number density of ions from the most 
probable value (in <7+i units), 

Ap(r) = A P ;(r) + Ap*_(r), (62) 

are related to the two OP's by 

g(k) = (Za(k) + 6(ife))0(k) + (Zb(k) - a(k))r}(k) (63) 
Ap(k) = (a(k) - b(k))4>(k) + (a(k) + 6(A:))^(k). 

The inverse relations are 

0(k) = ^ [(a(k) + b(k))q(k) + (a(k) - Zb(k))Ap(k)] (64) 

W = — J-j [(b(k) - ~a(k))m + (Za(k) + b(k))Ap(k)} • (65) 

On the low-temperature side of the structural line the dominant instantaneous states con- 
sist of the charge- and the number-density waves with the wavnumber k b and the amplitudes 
that satisfy the relation 

Ap(k b ) = Rq(k b ). (66) 

As seen from Eqs. (l63p - (l65p . the instability of Q MF associated with <j)(\t b ) (i.e. with fj(\<. b ) = 
0), corresponds to 

(Z~a(k b )+b(k b )Y 

For a(k b ) = b(k b ), i.e. in the RPM limit, R = 0, and (f)(k b ) oc q{k b ). When the instability of 
fl MF is induced by ^(k&) and the other OP vanishes, then 

R = (a(h) + b(k b )) 

(z~b(k b ) - d(k b )Y 

For Zb(k b ) = a(k b ) we have 1/R— 0, and fj(k b ) oc Ap(k fc ). 

For R <C 1 the charge- density waves dominate over the number-density waves in the 
dominant eigenmode. In such a case the regions of the excess positive charge are followed 
by regions of excess negative charge. The number densities in the charged regions are 
comparable. For R ^> 1 the number- density waves dominate, i.e. regions with excess 
number-density are followed by regions containing a smaller number of ions. Both the dense 
and the dilute regions are nearly charge-neutral for R ^> 1. 
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In real-space representation the OP fields are given by convolutions, 

0(r) = J dn[Ap;(n)o(r - n) - Ap* (n)6(r - n)], (69) 

with analogous expression for the field r](r), where a(r) and b(r) are inverse Fourier trans- 
forms of the functions defined in Eq. (|49p and (1501) . 

IV. CASE OF EQUAL SIZES 

A. MF spinodal and the wavevector of the most probable instantaneous state 

The phase separation (1531) occurs only for different sizes of the ionic species (5 7^ 0). 
In the size-symmetric case the disordered phase is unstable only with respect to periodic 
ordering in MF. Let us consider the periodic ordering, for which the MF spinodal is given 
by Eqs. (|56p and (|58p . For 5 = we obtain 

a aa = — + a+- (70) 

and 

n(A-s)(2 + s 2 ) 

a+ ~ = (l- a )4 • (71) 

The boundary of stability of the disordered phase and the associated wavevector reduce to 
the form 

24 cos k R 2 
t r( s ) = 72 s ' tanfctf = - — . (72) 

K R Kr 



The above form is identical to the one found for the RPM within the same approach[3|, 152 ]. 
The subscript R indicates that the MF spinodal and the wavenumber refer to the RPM. 
The charge asymmetry has no effect on the boundary of stability of the disordered phase on 



the MF level of our theory, as 



predictions of other theories 



17 



ong as the sizes are the same. This result agrees with MF 



59 



60 



60 



6l|. Beyond MF this property does not persist 



13, 



6l| due to the coupling between the two OP's in Q 



int- 



B. Eigenmodes 



In this case we obtain the following form of Eq. 



A=-* 9 nC,4k)— 2 ( Ts + ^r- ■ (73) 
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At the MF spinodal line (1721) A = 0. and in turn a(kn) = b{k R ) (see Eqs.gHJ) and (1501)). 
thus 

4>(k R ) = -L (Ap + (k R ) - Ap_(k fl )) (74) 

^(k fi ) = ^= (Ap + (k fi ) + Ap_(k*)) . (75) 

These relations are the same as in the RPM, although for Z / 1 0(k) differs from the 
charge- density amplitude fl6T|) . Beyond the spinodal line A^ 0, and the eigenmodes contain 
different proportions of Ap + and Ap_. 

V. THE MF SPINODAL ASSOCIATED WITH THE GAS-LIQUID PHASE SEP- 
ARATION 

A. Eigenmodes in the case of k — > 

Let us discuss the nature of the eigenmodes for k — > 0. From Eqs. (1491) and (150]) we find 
for arbitrary 5 and v 

a(0) = , Z , 6(0) = ; 1 . (76) 
In the long-wavelength limit the amplitudes of the eigenmodes assume the forms 

m = ^==f (Ap;(o) + w-(o)) (77) 

0(0) = ^L=(ZAp* + (0) - A/5* (0)). (78) 

Note that from Eqs. (1781) and flSTj) it follows that 0(0) oc g(0). The charge- neutrality condition 
((6]) is equivalent to 0(0) = 0. For 0(0) = we obtain 



v(P) = ^|^Ap(0). (79) 

In sec.IIIB we have shown that for — > the dominant eigenmode is f)(0). Eq. fl79l) shows 
that 7/(0) oc A/5(0). Hence, we find the usual phase separation into ion-dilute and ion-dense 
phases for all values of the asymmetry parameters, as expected. 
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FIG. 1: MF spinodal lines for the gas-liquid separation for 5 = 0.2 (A = 3/2). The solid and the 
dashed lines correspond to v > and v < respectively. From the top to the bottom, the solid 
lines correspond to v = 0.5, v = and v = 0.9. The upper and the lower dashed line correspond 
to v = —0.5 and v = —0.9 respectively. 

B. MF Spinodal lines 

In our MF approximation the spinodal line for the gas-liquid separation assumes a single 
maximum T*(( c ) at the critical point (cp) for any pair of 5 ^ and v. The MF spinodal lines 
are shown in Figs. 1,2 and 3 for 5 = 0.2, 5 = 0, 6 and 5 = 0.9 respectively, for several values 
of v. The spinodal lines assume the characteristic asymmetric shape as in the RPM l|, |2j]. 
The critical temperature and the volume fraction, T* and ( c respectively, are shown in Fig. 4 
as functions of the size asymmetry for several charge ratios. For a given charge ratio the 



dependence of ( c on the size asymmetry 



agrees 



qualitatively with the major trends found 



in Monte Carlo (MC) simulations 18|, [l9|, |20|, |22|, |44|. Namely, ( c is a convex function of 5, 
and its value is of the same order of magnitude as in Ref. 22|. In contrast, T* is a convex 
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2d, 



function of S for fixed u, whereas in MC a concave function was obtained 
For large asymmetries the value of T* is overestimated, whereas for 5 = we have T c 
in our MF approximation. On the other hand, the dependence of T* on v for fixed S agrees 



22]. 
= 



qualitatively with MC results 



22j . Namely, for 5 < the T* decreases with increasing charge 



asymmetry, and for 5 > the lines T*(S) corresponding to different charge ratios intersect 
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0.2 0.4 0.6 0.8 

FIG. 2: MF spinodal lines for the gas- liquid separation for 5 = 0.6 (A = 4). The solid and the 
dashed lines correspond to v > and v < respectively. From the top to the bottom, the solid 
lines correspond to v = 0.9, v = 0.5 and v = 0. The upper and the lower dashed line correspond 
to v = —0.5 and v = —0.8 respectively. 




FIG. 3: MF spinodal lines for the gasdiquid separation for (5 = 0.9 (A = 19). The solid and the 
dashed lines correspond to v > and v < respectively. From the top to the bottom, the solid 
lines correspond to v = 0.9, v = 0.5 and v = 0.0. The dashed line corresponds to v = —0.3. 



23 



each other (Fig.4 (upper panel) here and Fig. 9 in Ref. 22]). In the case of ( c , a correct 
dependence on the charge asymmetry is found for S > 0, whereas for 5 < our results do 
not agree with (partial) results obtained in Ref. 221]). 

In the case of 5 = 1 and v — > 1 (A = oo, Z — > oo), and for £ = 0(Z°), Eq. (1551) assumes 
the simple form 

T* = 3C(1-C) 2 - (80) 

The critical-point temperature, T* = 4/9, and density, ( c = 1/3, are both large compared to 
the corresponding values in the RPM, in contrast to the simulation results. Thus, some of 
the trends found in simulations are correctly predicted by the MF theory, while some other 
trends are not. Note, however that beyond MF coupling between the two OP's may have a 
significant effect on the spinodal. In particular, beyond MF T* > is found for 8 — [3], 0] 
as a result of the coupling between the two OP's in Q int (see Eq. (l25]) ). Significant effect of 
the coupling between the two OP's on the spinodal should be expected also for 5 7^ and 
v 7^ 0. Indeed, a qualitatively correct dependence of the location of the cp on v is found 
for 5 = in the CV approach jl7j ]. beyond MF. As mentioned in sec.IIA, the CV approach 



and our mesoscopic theory are closely related [49] . Correct trends are also found when ion 
pairing is explicitly taken into account Q • 



VI. THE MF SPINODAL ASSOCIATED WITH PERIODIC ORDERING OF IONS 

The structural line (MF spinodal) associated with the periodic ordering of the instan- 
taneous states, as well as the most probable periodic structures depend qualitatively on 5 
and v. The key role is played by the size asymmetry. Three regimes of the size asymmetry 
can be distinguished, although there are no sharp boundaries, but rather smooth crossovers 
between them. The approximate ranges of the first two regimes are 8 < 0.4 ( A < 2.3), and 
0.4 < 5 < 0.9 (2.3 < A < 19). The qualitative dependence on v is found in the second regime 
only. Third regime corresponds to the asymptotic case of 5, v — > 1, i.e. to extreme charge- 
and size asymmetry. In the other (unphysical) extreme, 8, v — > —1, the periodic ordering is 
suppressed. Detailed description of each case is given in the following subsections. 
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FIG. 4: Critical temperature (top) and critical volume fraction (bottom) as functions of the size 
asymmetry 5, for three values of the charge asymmetry v. Solid, dashes and dash-dotted lines 
correspond to v = 0.9, 0.5 and respectively. 

A. Small size asymmetry 

1. MF spinodal lines and wavevectors in the dominant instantaneous states 

The MF spinodal in the case of the small size-asymmetry 5 < 0.4 (A < 2.3) is qualitatively 
the same as in the RPM, i.e. along the spinodal T* is a monotonically increasing function 
of ( (Fig. 5). Somewhat surprising is the fact that for larger charge-asymmetry the ordering 
occurs at higher temperatures. It is not clear whether this tendency persists beyond MF. 
The wavenumber characterizing the period of the oscillations of the OP field in the ordered 
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FIG. 5: MF spinodal lines for the transition to the ordered phase for 5 = 0.2. From the top to the 
bottom lines v = 0.9, v = 0.7, v = 0.5 v = 0.0, v = —0.5, v = —0.7 and v = —0.9. Temperature 
T* and the volume fraction of ions £ are in dimensionless reduced units defined in Eqs. (|ll|) and 
(fl"2"j) respectively. 



phase is k b > ir/2, and its value is almost independent of £ (see Fig. 6). The size of the regions 
of increased or depleted densities of the two ionic species, ~ 7r/fct, < + is comparable 
to the size of the ions, and the ordered phase should be identified with a standard 'hard' 
ionic crystal. According to the discussion in sec. III. B, the relevant OP field is (p{k^) for the 
whole temperature range. The amplitudes of the number- and the charge-density waves in 
the dominant instantaneous state are related according to Eqs.( l66|) - (!67|) . The amplitude 
ratios R for S = 0.2 are shown in Fig. 7. 

The above results indicate that the behavior of the PM for S < 0.4 is qualitatively the 
same as in the RPM. By analogy with the RPM, we expect that the transition to the ionic 
crystal is fluctuation-induced first-order, and occurs at significantly lower T* and higher £ 



compared to the structural line [37]. 
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FIG. 6: The wavenumber fcj corresponding to the ordering of ions along the MF spinodal lines 
shown in Fig. 5 for 5 = 0.2. From the bottom to the top lines v = 0.9, v = 0.7, v = 0.5 v = 0.0, 
v = —0.5, v = —0.7 and v = —0.9. kj, is in cr^i units and £ is the volume fraction of ions. 

2. Eigenmodes for very small asymmetries 



Let us determine the behavior of the eigenmodes in the RPM limit. For 5, v — > we find 
from Eq.flMD that A^0, and 



1 



where 



a(k) ~s*-+o ^(1 + y) + 0(5 2 , Sv, v 2 ) 
b(k) ~ s ^^L(i-y) + 0{8 2 ,8v,v 2 ) 



The eigenmodes (j47j) and (j48j) take the forms 



0(k) 



v/2 



Ap + (k) - Ap_(k) + y(Ap + (k) + Ap_(k) 
1 

"72 



1 - i/)g(k) + (y - i/)Ap(k) + 0(<5 2 ,<5z/,z/ 2 ) 



(81) 



(82) 



(83) 
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FIG. 7: The ratio R in the relation (|66p between the amplitudes of the number- and charge-density 
waves in the dominant eigenmode, along the MF spinodal lines shown in Fig. 5 for 5 = 0.2. From 
the top to the bottom lines (on the right) v = 0.9, v = 0.5, v = 0.0, v = —0.5 and v = —0.9. 
Both R and the volume fraction of ions £ are dimensionless. Note that \R\ < 1, indicating periodic 

ordering into positively and negatively charged regions of the spacial extent 7r/fc& (in er_| units), 

where kb is shown in Fig. 6. 



and 



i 



V2 



Ap + (k) + Ap_(k)-y(Ap + (k)-Ap_(k)) +0(5\5^u 2 ) (84) 

= 4|[Ap(k) -yg(k)] +0(5 2 ,^ 2 ). 

For infinitesimal asymmetry parameters, Ap and g yield infinitesimal contribution to <fi and 
to 77 respectively. The OZ equations are thus decoupled for the eigenmodes <fi and rj that 
very weakly deviate from the charge- and the number densities respectively. This result 
agrees with the predictions of Ref . pj . 
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B. Moderate and large size asymmetry 

For 5 > 0.4 (A > 2.3) the structural line associated with the periodic ordering of the 
instantaneous states assumes a single maximum for ( « 0.15. Unlike the case of 5 < 0.4, 
at higher volume fractions the periodic ordering is less favorable. For the nonmonotonic 
structural lines two cases can be distinguished, according to a different behavior of the 
wavenumber of the dominant instantaneous states kb- 



1. k b > tt/2 

Such large wavenumbers are found for the moderate size asymmetry, 0.4 < 5 < 0.67, in 
the case of 5u > (larger charge at the larger ion). For 5 = 0.6 and v > the structural 
lines and the corresponding wavenumbers kb are shown in Figs. 8 and 9 respectively. Except 
for very small (, the wavenumber is kb > tt/2, and increases with increasing (. Increasing 
kb indicates decreasing period of oscillations of p Q (x), in agreement with expected smaller 
nearest-neighbor distances for larger density. For kb > tt/2 the dominant field is (sec.IIIB). 
The behavior of R, shown in Fig. 10, is similar as in the case 5 < 0.4. 



The analysis of the cluster formation in Ref.|46j was performed for this range of param- 
eters. Namely, the case of 5 = 0.639, and v = 0.85 (monovalent), v = 0.7 (divalent) and 
v = 0.6 (trivalent counterions) was examined in detail. The volume fractions and temper- 
atures were (C,T*) = (0.0185, 0.157), (0.0174, 0.079) and (0.017,0.0525) for mono-, di- and 
trivalent counterions respectively. All cases correspond to the low-T* side of the structural 
line, where the periodic instantaneous structures dominate. Only for trivalent counterions 
kb < vr/2 and the relevant OP is rj, but this fact is not sufficient to explain the differences 
found in the short-range ordering into clusters. Further studies beyond the stability analysis 
of Q MF are required to verify whether the coarse-grained description can explain the cluster 
formation. 



2. k b < it/2 

kb < tt/2 is found in the case of 0.4 < 5 < 0.67 for v < (larger charge at the smaller 
ion), as shown in Fig. 11. 
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FIG. 8: MF spinodal lines of the transition to the ordered phase for 5 = 0.6 and v > 0. From the 
top to the bottom lines v = 0.9, v = 0.7, v = 0.5 and v = 0.0. Temperature T* and the volume 
fraction of ions £ are in dimensionless reduced units defined in Eqs. ([ll|) and (|12|) respectively. 
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FIG. 9: The wavenumber fc& corresponding to the ordering of ions along the MF spinodal lines 
shown in Fig. 8 (for 5 = 0.6 and v > 0). From the bottom to the top lines v = 0.9, v = 0.7, v = 0.5 
and v = 0.0. fc& is in cr^ 1 : units and £ is the volume fraction of ions. The dashed and dash-dotted 
lines correspond to kb = tt/2 and ir respectively. Note that fc& < 7r/2 for very small values of £. 
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FIG. 10: The ratio R (|66p between the amplitudes of the number- and charge-density waves in the 
dominant eigenmode for 5 = 0.6 and v > 0. The solid, dashed, and dash-dotted lines correspond to 
v = 0.9, v = 0.5 and v = 0.0 respectively. Note that \R\ < 1, and periodic ordering into positively 
and negatively charged regions of a size ir/kb, with kb shown in Fig. 9, takes place. Both R and the 
volume fraction of ions £ are dimensionless. 
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FIG. 11: The wavenumber kb corresponding to the ordering of ions along the MF spinodal lines 
shown in Fig. 12 (for S = 0.6 and v < 0). From the top to the bottom lines v = —0.5, v = —0.7 and 
v = —0.8. kb is in units and £ is the volume fraction of ions. At the dashed line kb = vr/2. 
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0.2 0.4 0.6 0.8 

FIG. 12: Solid lines represent the MF spinodal lines of the transition to the ordered phase for 
5 = 0.6 and v < 0. The dotted lines are the MF spinodals corresponding to the gas-liquid phase 
separation associated with kb = 0, discussed in sec.V. From the top to the bottom lines v = —0.5, 
v = —0.7 and v = —0.8. Temperature T* and the volume fraction of ions Q are in dimensionless 
reduced units defined in Eqs. lfTTj) and ([12]) respectively. 

For 5 > 0.67 (A > 5) the wavenumber of the dominant eigenmode of C®Jkb) is kf, < tt/2 
for all values of v (Fig. 13). The extent of the phase-space region corresponding to periodic 
ordering of the instantaneous states is similar to that of the case of 0.4 < 5 < 0.67 (Figs. 14 
and 15). 

When the difference between v and 5 is not too large, the kj, initially increases with 
increasing ( (Fig. 11 (upper curve) and 13 (two upper curves)). When fc 6 increases, the 
structural line and the MF spinodal of the gas-liquid separations are well separated. The 
amplitude ratio is R = 0(1) and varies slowly with (. At a value of ( depending on 8 
and u, kf, starts to decrease, and at some point a rapid decrease to k\, = occurs. At the 
corresponding value of ( = (l (Lifshitz point) the MF spinodals merge together and become 
identical for C > Cl (Figs. 12, 14 and 15). The amplitude ratio R ( |66l) increases rapidly, 
when a rapid decrease to kb = occurs. For a large difference between the charge and the 
size asymmetry, the kb decreases for the whole range of ( (see the two lower curves in Fig. 11 
and the three lower curves in Fig. 13). 
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FIG. 13: The wavenumber kb corresponding to the ordering of ions along the MF spinodal lines 
for 5 = 0.9. From the top to the bottom lines v = 0.9, v = 0.7, v = 0.5, v = 0.0 and v = —03. 
kf, is in units and £ is the volume fraction of ions. The horizontal dashed line corresponds to 
k b = it/2. 
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FIG. 14: MF spinodal lines (solid) for the transition to the ordered phase for 5 = 0.9 and v > 0. 
Dotted line represents the MF spinodal line for the gasdiquid separation. From the top to the 
bottom lines v = 0.9, u = 0.7 and v = 0.5. Temperature T* and the volume fraction of ions £ are 
in dimensionless reduced units defined in Eqs. (|lip and (|12p respectively. 
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FIG. 15: MF spinodal lines (solid) for the transition to the ordered phase for 5 = 0.9 and v < 0. 
Dotted lines are the MF spinodal lines for the gas-liquid separation. From the top to the bottom 
lines v = 0.0 and v = —0.3. Temperature T* and the volume fraction of ions £ are in dimensionless 
reduced units defined in Eqs. lfTTj) and (fT2|) respectively. 

For R ^> 1 the amplitude of the charge-density is much smaller than the amplitude of 
the number-density wave (sec. HID). Thus, the ordering into alternating oppositely charged 
regions is suppressed. For v < a large charge at the small ion is compensated by a large 
number of large ions having small charges; therefore the unit cell of the ordered structure is 
large, and may contain voids. In crystals containing voids as parts of the structure the size 
of the unit cell is significantly larger than the sum of radii of the two ionic species. Such 
voids should be distinguished from vacancies resulting from thermal fluctuations, which are 
present in the ionic crystals at low densities. 

In this case the favorable ordered structures are much more complex than in both, the 
RPM-like systems and the colloid-like systems. Whether such complex structures may cor- 
respond to stable phases in the PM with moderate size- and charge asymmetry remains an 
open question, because the binodal associated with the periodic ordering may be preempted 
by the binodal associated with the phase separation. 

For a large difference between the charge and the size asymmetry, the MF instability with 
respect to ordering of ions in periodic structures occurs only for small values of ( and T* 



34 




t i I i I i I i I i I i I d 
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FIG. 16: The ratio R (|66p between the amplitudes of the number- and charge-density waves in the 
dominant instantaneous state for 5 = 0.6 and v < 0. The upper and the lower curve correspond to 
v = —0.8 and v = —0.5. Note that R > 1, and regions of increased and depleted number density 
are formed. The size of these regions is n/kb, with kb shown in Fig. 11. Both R and the volume 
fraction of ions £ are dimensionless. In the case of the upper curve both the dense and the dilute 
regions in the dominant eigenmode are nearly charge neutral for £ > 0.1. 

(Figs. 12, 14, and 15). For v5 — > —1 the ordering is suppressed entirely. An infinite number 
of large ions is required to neutralize the infinite charge at the small ion for v — — 1 . Clearly, 
formation of an ordered structure with a finite period is not possible in this (unphysical) 
case. 

Finally, we should note that at large volume fractions the effects of ordering of hard 
spheres become important, and such effects are not accounted for by this theory. 

C. Case of very large asymmetry, Su — ► 1 

The maximum of the MF spinodal T*(Q increases when v w 8 — > 1, and for v ~ 5 > 
0.9 the maximum increases very rapidly. This increase of the phase-space region where 
a majority of the instantaneous states are the periodically ordered can be observed by 
comparing Figs. 13, 18 and 19. In the latter plot the spinodal line and the wavevector are 
shown for Z = 10 4 (y = 0.9998) and A = 10 3 (5 = 0.998), corresponding to highly charged 
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FIG. 17: The ratio R (|66[) between the amplitudes of the number- and charge-density deviations 
from the average values in the dominant instantaneous state for 5 = 0.9. From the bottom to the 
top lines v = 0.9, 0.5, 0, and —0.3. At the dashed line R = 1. Both R and the volume fraction of 
ions C are dimensionless. 



colloids with counterions of a microscopic size, considered in many experimental studies. 

The case of extreme asymmetry was studied in Ref.Jsj] for s = O(Z ) = 0(Z~ r )), 
i.e. for infinite dilution of colloid particles. In the limiting case the spinodal assumes the 
asymptotic form 8| 



3 cos(2fe 



c 



k c 



Z( 



where 



tan(2/c 



ci 



1 



(85) 



(86) 



The subscript C indicates that the MF spinodal and the corresponding wavevector corre- 
spond to the colloid limit, 5, v — ► 1, with the infinite dilution of the colloid particles for 
Z — > oo. 

Note that kc = kn/2 < n/2, and according to the discussion in sec.IIIB the eigenmode (j> 
is relevant on the high-T* side, whereas the eigenmode r] is relevant on the low-T* side of the 
line ( 1521) . The asymptotic analysis in Ref.js| indicates that the phase behavior is determined 
only by the field Ap + . Here we shall verify if the dominant eigenmode indeed reduces to 
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FIG. 18: The MF spinodal lines (solid) and the corresponding asymptotic behavior for v = 5 — > 1 
(dashed). The latter has the form T* = 12.3C/(1 - S 4 ) (Eq.38 in Ref. 8] and Eq.flU]) here). From 
the bottom to the top lines 5 = v = 0.95, 0.98, 0.996 and 0.9995 (A = Z = 39, 99, 499 and 3999) 
respectively. In the last case the solid and the dashed lines are indistinguishable for the given range 
of T*. 

Ap + . For Z — > oo and for C ++ {k) ^ we obtain the approximate form of Eq. (j4"5"]h 

2C (k) 2 

B = \C ++ (k)\ - sign(C ++ (k))C^(k) + + ~ / + 0(Z~ 2 ). (87) 

In the limit 8 = v = 1 the line ((52J assumes the form T* = ^p^. From (JU]), §7$ and the 
above we obtain the asymptotic behavior of a and b for Z —>■ oo, 

1 ^ °(J0 if H?^ < T * < T c(C) or 7£(C) < T* < 3cosfc 



3(Jfc) = <~ fc2 CVW *cvs,^ - fc2 (g8) 

O ( J otherwise , 



with = O J and = 1 — 0(-^J in the first and in the second regime respectively. 
The eigenmodes obtained from Eqs.(j47"l) and ( )48l) assume the forms 



o(k . Ap + (k) if ^<T*<T£(C) or T5(C)<T*<^|A, 

Ap_(k) otherwise, 
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FIG. 19: The MF spinodal line (top) and the corresponding wavevector of the dominant eigenmode 
(bottom) for 5 = 0.998 and v = 0.9998. The dashed line is the asymptotic result of Ref. S\. For 
C <C 1, kf, ~ 1.23 in agreement with the asymptotic result of Ref.[8j]. T* is in the dimensionless 
units (Eq, (fTT|) ). and kb is in units of a^L Rt 2/a + . 



and 



Ap_(k) if 



3 cos k 
k 2 



< T* < T£(C) or T*(C) <T* < 



3 cos k 
k 2 ' 



(90) 



Ap + (k) otherwise, 

where the neglected contributions are 0(1/Z). The dominant eigenmodes for different parts 
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FIG. 20: Dominant eigenmodes for 8, v — * 1 (A, Z — * oo) in different parts of the phase diagram. 
As discussed in sec.IIIB, the fields rj and eft dominate below and above the dashed line respectively. 
The relation between the eigenmodes and Ap a is given in Eqs. (|89p and (|90p . The neglected 
contributions to the OP's are 0(1 /Z). 

of the ((,T*) phase diagram, obtained from Eqs. (1891) and ( 1901) . are shown in Fig.20. Below 
the MF spinodal line, i.e. for T* < Tq((), each dominant OP reduces to Ap + (kc)- The 
asymptotic analysis of Ref. 8| is thus fully consistent with the present, more complete 
approach. 



VII. SUMMARY 



We performed the MF stability analysis of the disordered phase in the PM for the whole 
range of size and charge asymmetry of the ionic species. At the true instability the average 
distribution of ions undergos ordering. As we argue in sec.IID, the instabilities of the 
disordered phase found in MF in fact represent structural lines. The structural lines are 
associated with ordering of the majority of the instantaneous states, i.e. with pretransitional 
effects. Our results show that periodic ordering of the instantaneous states depends mainly 
on the size asymmetry 5. Qualitative dependence on the charge asymmetry v is found only 
when the size asymmetry is sufficiently large. 

The strongest tendency for periodic ordering is found for 8, v — > 1; in this case the 
pretransitional effects are present for a very large temperature range (Fig. 19, top). The 
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period of p a (r), 27i/kb ~ 2cr + , decreases with increasing ( (Fig. 19, bottom). More detailed 
analysis in Ref. [8J shows that for small £ the most probable structure has the form of the 
bcc crystal formed by the large ions that are immersed in the cloud of the counterions. 
When both asymmetry parameters decrease the temperature range for which the periodic 
ordering occurs decreases very rapidly (Figs. 18, 14, 8). For 8, v — 0.9, T* at the maximum 
of the structural line is 500 times lower than for 5 = 0.998, v = 0.9998 (Figs. 14 and 
19). Physical systems corresponding to this range of asymmetry parameters include highly 
charged colloids in the presence of one kind of counterions and no coions. 

In the other extreme case of bv — > — 1 the periodic ordering is suppressed entirely. This 
case is unphysical, and will not be discussed here. 

For 0.9 > 5 > 0.4 the periodic ordering depends on the size asymmetry rather weakly, 
as long as v is sufficiently large (Figs. 8 and 14). However, when v decreases and becomes 
negative (the charge at the smaller ion is not sufficiently small), T* and ( corresponding to 
the ordering both decrease to very small values (Figs. 12 and 15). Nearly neutral clusters 
followed by voids are formed in this case, and the period of the structure, 27r /kb, is an 
increasing function of £ (Figs. 11, 13). Moreover, the corresponding structural line and the 
MF gas-liquid spinodal are close to each other on the phase diagram (Figs. 12 and 15). In 
future studies the binodal lines should be determined in MF and beyond to verify whether 
the complex ordered structures with large periods correspond to stable crystal phases, or 
whether they reflect a tendency for formation of aggregates in the fluid phase. Physical 
systems corresponding to this range of asymmetry may include globular proteins, organic 
ions, and some room-temperature ionic liquids. 

For the large and moderate size asymmetry discussed above (S > 0.4) the structural line 
assumes a pronounced maximum for <^ ^ 0.15 (Figs. 8, 12, 14, 15, 18). It is remarkable that 
the periodic ordering induced by the Coulomb interactions is most efficient for small volume 
fractions, whose range, 0.1 < ( < 0.3, is almost independent of 5 (Figs. 8, 12, 14, 15, 18). 

When 5 further decreases, the maximum of the structural line becomes less and less 
pronounced, and finally dissappears for 5 ~ 0.4 (A < 2.3). T* at the structural line becomes 
a monotonically increasing function of £, and the periodic ordering is more efficient at 
large volume fractions (Fig. 5). The period of the ordered structure is nearly independent 
of ( (sec. VIA). The same qualitative behavior is found in our MF approximation for the 
RPM 52|. By analogy with the RPM 33] we expect that beyond MF the ordered phase 
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corresponds to a hard ionic crystal. The crystallization line is expected at much lower T 



371 ] . Physical systems for this range of 



and much higher £ compared to the structural line 
asymmetry include molten salts and electrolytes. 

Finally let us focus on the separation between two uniform phases. We found that some 
of the effects of the size- and charge asymmetry are correctly predicted already on the MF 
level of our theory (sec.V). However, on the MF level the gas-liquid separation is preempted 
by the instability with respect to the periodic ordering, except from the case where such 
ordering is suppressed by geometrical constrains. This fact agrees with earlier MF results 



for the RPM 
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52 



571 ]. In fact, on the MF level no phase separation is found in the RPM; 

58] . Studies 
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the phase separation is induced by the charge-charge correlations |41| . 
beyond MF are required to determine for what parts of the parameter space ((, T*; (5, i/) 
the phase separation is indeed preempted by the crystallization, and how the critical point 
varies with 5 and v. Returning to the clustering associated with the phase separation we 
should note that phenomena occurring at the very short length-scale cannot be accurately 
described by the coarse-grained theory. It is not clear to what extent the clustering of ions is 
associated with the short-distance behavior of the correlation functions, and to what extent 
it is induced by collective phenomena. Studies beyond the stability analysis of Q MF may 
shed light on this question. 

Short-range interactions of different origin (including solvet-induced effective interactions 
between solute molecules) also play an important role for the phase behavior [3J. Such 
short-range interactions can be included in the field-theoretic approach by supplementing 
the energy contribution to fl MF in Eq.flZJ) with additional terms. 
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VIII. APPENDIX A 



The quantities a a p in Eq. (!20|) are obtained from the explicit expressions for the chemical 



potentials in the hard-sphere mixture derived in Ref. 
the following expressions 



48| . Direct differentiations lead to 



a_H 



1 ixr\ 
Po+ 6 



.i-C ' (i-C) 2 

with the analogous expression for a , 



^ 6r 2 + Y + 15r+X + r% s + 18r 2 + X 2 + 6r 3 + XY + 9riX 3 



(i-C) ; 



(91) 



8 2r+r_(6 + r + r_)X + 8r 2 + r 2 _Y 18r 2 + r 2 _ X 2 + Qr 3 + r'iXY 9r 3 + r^X 3 

1 ~ ~ 1 ~ ~ 1 — ~ m - (92) 



.i-C ' (i-C) 2 (i-C) 3 (i-C) 4 

where s, ( and r ± are defined in Eqs. (|T3j) . (|T2|) and (124]) . and following Ref. 
introduced the notation 



X 



71 

6' 



S Pl+ + rlpl-) = *(l + S' - 25v) 
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Y =Q( r +Po + + r -Po-) = s ( 1 -^)■ 



we have 

(93) 
(94) 



IX. APPENDIX B 



The functions b and d introduced in Eq. ( IBTl) are defined by 



b = (a ++ Z 1 + a Z) cos k cos(k5) + (a ++ Z l — a Z) sinfcsin(M) + 2a + _ cos/c, (95) 



1 (1 + 2C) 2 



d = a i + a — a, = ; — > 0, (96) 

++ + " Po + Po- (1 - C) 4 1 ^ 
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